Multireference Correlation in Long Molecules with the Quadratic Scaling Density 

Matrix Renormalization Group 
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We have devised and implemented a local ab initio Density Matrix Renormalization Group 
(DMRG) algorithm to describe multireference correlations in large systems. For long molecules 
that are extended in one of their spatial dimensions, we can obtain an exact characterisation of 
correlation, in the given basis, with a cost that scales only quadratically with the size of the system. 
The reduced scaling is achieved solely through integral screening and without the construction of 
correlation domains. We demonstrate the scaling, convergence, and robustness of the algorithm in 
polyenes and hydrogen chains. We converge to exact correlation energies (with 1-lOpEh precision) 
in all cases and correlate up to 100 electrons in 100 active orbitals. We further use our algorithm to 
obtain exact energies for the metal-insulator transition in hydrogen chains and compare and contrast 
our results with those from conventional quantum chemical methods. 



I. INTRODUCTION 

The electronic structure of a chemical system features 
two types of electron correlation. The first is nondynamic 
correlation. This is associated with the correlation of 
electrons in nearly degenerate valence orbitals. The cor- 
rect description of nondynamic correlation is necessary 
to establish the qualitative features of chemical bonding. 
The second is dynamic correlation. This is associated 
with excitations from valence degrees of freedom into the 
many non-bonding orbitals. The multiple weak excita- 
tions are responsible for establishing the detailed, quan- 
titative structure of the electronic wavefunction. 

In general, a correct description of strong nondynamic 
correlation in large systems is very difficult to obtain. 
When nondynamic correlation is important (e.g. during 
bond breaking) , a single determinant or electronic config- 
uration does not provide the correct qualitative structure 
of the wavefunction. Instead, the delicate balance in the 
valence degrees of freedom between the kinetic energy, 
which favours delocalisation, and the Coulomb energy, 
which favours localisation, results in competing electronic 
configurations and the correct electronic structure con- 
tains contributions from multiple determinants with sig- 
nificant weights. Complete- Active-Space Self- Consistent- 
Field (CASSCF) theories correctly describe this type 
of structure by expanding the wavefunction in the com- 
plete space of the optimised valence (or "active" ) degrees 
of freedom, but do so at the cost of a factorial scaling 
with the number of active electrons. Such calculations 
with more than O(10) electrons remain extremely diffi- 
cult at this time. Despite the impressive progress in lo- 
cal Generalised Valence Bond and Coupled Cluster (CC) 
Theories (e.g. iiSIillillElElElElEI) 
which provide some capacity to break e.g. single bonds, 
such approaches do not possess the flexibility of a true 
multireference theory. The fundamental challenge there- 
fore remains to find a multireference electronic struc- 
ture method that is sufficiently flexible to correctly de- 
scribe nondynamic correlation, yet which exhibits a non- 



factorial scaling, and can thus be applied to large sys- 
tems. 

In this work, we will adopt the more modest goal 
of answering the question of how to describe nondy- 
namic correlation in systems which are large in only 
one out of their three spatial extents. In quasi-one- 
dimensional systems the physics that is familiar from 
three-dimensions, is notably modified. This is illustrated 
by the organic electronic materials (e.g. conjugated or- 
ganic polymers and carbon nanotubes) which exhibit un- 
usual interacting electron effects, arising from coupled 
quasi-one-dimensional motions of many electrons along 
the conjugated 7r-backbone. As a simple example, in lin- 
ear polyenes, electron-electron interactions contribute to 
make the lowest excited singlet state (the 1 2A g state) 
one of double-excitation nature, rather than the singly 
excited HOMO — > LUMO state as one would expect in 
a single particle picture 0, fl6l Il7| . A more extreme ex- 
ample of this occurs in metallic nanotubes at low temper- 
atures, where at low energies there are no singly excited 
states; all low-energy excitations are of collective nature, 
and the qualitative electronic structure is of Luttingcr 
liquid form HI E1I23 . 

Here we demonstrate that the Density Matrix Renor- 
malization Group (DMRG), which we and others have 
recently been inv estigating in qu antum chemistry [2lJ, 123, 

11 pi pi El El ElM IM .Hi 111 111 Hip EE IH HI 

133. l40l lilj , provides a solution to the question of how to 
flexibly and efficiently describe nondynamic correlation 
in systems that are large in one of their three spatial di- 
mensions. Our analysis shows that the DMRG behaves 
as a local, multireference, size-consistent/size-extensive, 
and variational theory. From the intrinsic locality of the 
DMRG ansatz we formulate a DMRG algorithm, denoted 
for convenience as LDMRG, that scales only quadratically 
with the size of the system, without any need for an ar- 
tifical imposition of orbital domains. The multireference 
nature of the ansatz also eliminates any need for sepa- 
rately localised occupied and virtual orbitals. Using this 
algorithm, we carry out numerically exact DMRG calcu- 
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lations for long molecules, including polyenes in the Tr- 
active space and metallic and insulating hydrogen chains 
where we correlate up to 100 active electrons in 100 active 
orbitals. 

The structure of our discussion is as follows. We first 
introduce the DMRG wavefunction ansatz in section ITT1 
There we discuss its multireference, size-consistent and 
size-extensive, variational, and local properties, and the 
implications for the design of a local DMRG algorithm. 
In section ITTT1 we show how a simple screening of integral 
amplitudes results in a robust and naturally quadratic- 
scaling DMRG algorithm. In section llVl we present cal- 
culations on hydrogen molecular chains and polyenes in 
the 7r-active space and demonstrate the size-extensivity, 
computational scaling, and convergence of the LDMRG 
algorithm. As a difficult test of nondynamic correlation, 
we further carry out calculations on the metal-insulator 
transition in hydrogen chains for both symmetric and 
asymmetric bond stretching, and compare our results 
against existing quantum chemical methods (section Ej. 
Finally, our conclusions are presented in section IVll 



II. THE DMRG ANSATZ 

A. DMRG and Matrix Product States 

In previous work [29|. we have described the Renor- 
malization Group formulation of the DMRG al gori thm. 
However, as related by Ostlund and Rommer |42l |43| . 
and subsequently developed by other authors (see e.g. 
0113 ), the DMRG is also fruitfully analysed from the 
viewpoint of the underlying wavefunction ansatz, the Ma- 
trix Product State (MPS). Such a formulation will be 
convenient for our present discussion and we will recall 
the main points below; for a full presentation, we refer 
the reader to the excellent review by Schollwock |4q . 

Consider an iV-particle system in a state |\&) spanned 
by k orbitals. In occupation number representation, 
can be expanded as 

l*)= E ^nin a ...n>l"3 •••»*) ( IL1 ) 

nw2—nk 

where rii — 0,1; Y^i n i = N. We now decompose the 
high dimensional coefficient tensor tp into a chained ma- 
trix product via repeated singular value decompositions 
(SVD). For example, if there are only two orbitals, a sin- 
gular value decomposition yields 

^n^E^r^r (n.2) 

i 

where R™ 1 and R™ 2 are the singular vectors and <x are the 
singular values. Similarly, for three orbitals, ip nin2n3 can 



be decomposed via two singular value decompositions as 

i 

= j2 Rr r^R%R? (n.3) 

ij 

where in the SVD of Si, all singular values have modu- 
lus one, since Si is an orthogonal matrix. In this way, 
through repeated SVDs, the fc-dimensional coefficient 
tensor can be decomposed as a chain of matrix products, 

Vw.« fc =Tr{R»iR" 2 ...<T...R"*}. (II.4) 

So far the decomposition is exact since the R matrices 
arc full rank and will grow increasing large as the number 
of orbitals grows. 

The Matrix Product State which underlies the DMRG 
algorithm, arises by truncating the maximum dimension 
of the R matrices to be at most M x M, and thus with 
this restriction, we write the MPS as 

|#)= E Tr {R ni R" 2 . . .cr . . . R" fc } |mri2 ■ • ■ n k ). 

n 1 n 2 ...n k 

(II.5) 

We now establish the relationship between the MPS 
and the usual formulation of the DMRG algorithm. Re- 
call that any point in a DMRG sweep, the orbitals are 
partitioned into two blocks: a left block (spanning or- 
bitals 1, . . . , / say) and a right block (spanning orbitals 
g = f + 1, . . . , k). Through successive renormalisation 
transformations we obtain an adaptive many-body basis 
of dimension M to span the orbitals 1, . . . , /; let us de- 
note these many-body states by |Z/). First we enlarge 
the left block by adding the next orbital to give a super- 
block with an associated space {|//)} ® {|^)}- Next we 
renormalise this space to form a new many-body basis 
for the enlarged block spanning orbitals l,...,g, 

as 

l^=E</IW (H.6) 

where the rows of the matrix R™ 9 are the M eigenvectors 
of the density matrix of the superblock 1, .. .,g. After 
successive renormalisations, we see that the renormalised 
states take on a matrix product form, e.g. 

h) = E^iw 

gnu 

= E R l h g R 7f\ l f U 9 n h) 
fgn g n h 

{11.7) 

where each R Ki matrix is truncated to have maximum 
dimension M x M, 

To complete the identification of the underlying 
DMRG wavefunction with the Matrix Product State, 
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we introduce the corresponding renormalised many-body 
states \r g ) which span the orbitals g = f + 1, . . . , k. In 
the tensor-product space of the left and right blocks, we 
can write the full wavefunction in the form 

w = Evjw- (n.8) 

Performing an SVD, we obtain 

i*> = Ei ; ~>/i f />- ( IL9 ) 
/ 

Substituting in the matrix product decomposition of the 
DMRG many-body basis for the left and right block 
basis states from (|II.7ll in eqn. (|II.9(1 we identify the 
DMRG wavefunction with the Matrix Product State 
(fTl5)l . Consequently, the DMRG can be viewed as a self- 
consistent optimisation algorithm for the Matrix Prod- 
uct State where the renormalisation matrices R Mi which 
parametrise the ansatz are determined one by one from 
the density matrices of the blocks after each blocking step 
in a DMRG sweep. The number of retained states in the 
DMRG M thus coincides with the dimensionality of the 
matrices that parametrise the MPS. We note that the 
position of er in the Matrix Product State corresponds 
to the point of division between left and right blocks in 
the DMRG algorithm. In principle, the DMRG wave- 
function varies with different block partitionings along a 
sweep, but in practice, the variation is quite small. 



B. DMRG as a local, multireference, variational, 
size-consistent ansatz for long systems 

Starting from the perspective above, let us summarise 
some features of the DMRG/MPS ansatz. 

1. Variational: Since we can associate a wavefunc- 
tion with any DMRG block configuration, and a 
DMRG energy is evaluated as an expectation value 
of such a wavefunction, the energies appearing in 
the DMRG procedure are strictly variational. 

2. Multireference: It is clear that the Hartree-Fock 
reference has no special significance in the DMRG 
state, and in particular, we do not order or rank 
excitations relative to a single reference state. Fur- 
thermore, in contrast to selected Configuration In- 
teraction (CI) theories, none of the coefficients of 
expansion ip nin2 ... nk are restricted to be zero. 

3. Size-consistency: Within a physical ordering of the 
orbitals on the DMRG lattice, the Matrix Prod- 
uct State for two widely separated systems fac- 
torises into the product of Matrix Product States 
for each system separately. To see this, first ar- 
range the orbitals into left and right blocks, with 
the left block containing orbitals of the first system, 



and the right block containing orbitals of the sec- 
ond system. Since there is no coupling, the Matrix 
Product State for the total block configuration is a 
product |$) = \l)\r), where |Z) is a Matrix Product 
State for the first system considered alone (without 
changing the orbital ordering) and similarly for \r). 
Consequently, the DMRG energy is size-consistent. 

4. Locality and compactness: The number of varia- 
tional parameters in the Matrix Product State is 
0(M 2 k) and its correlation length is determined by 
M. Thus in any system with a finite quantum (i.e. 
off-diagonal) correlation length along the DMRG 
lattice, we can obtain a given accuracy in the en- 
ergy per unit site with constant M, independent of 
the size of the system. In such cases, for a given ac- 
curacy, the number of variational parameters in the 
DMRG scales linearly with the size of the molecule. 
The restriction to given M determines the finite 
correlation length that is captured by the ansatz; 
there is no need to a priori impose any orbital do- 
mains. Thus the DMRG is a naturally local scaling 
ansatz, and so long as the determination of the en- 
ergy is also performed with an account of locality 
(e.g. through screening or multipole expansion) a 
low-order scaling correlation theory arises. Indeed 
this is the basis of the quadratic scaling algorithm 
in the next section. Note that a finite correlation 
length implies only that we are away from a quan- 
tum critical point; such wavefunctions need not be 
close to the Hartree-Fock reference in any sense, 
as is indeed the case for systems with strong inter- 
actions. Thus, the local correlation nature of the 
DMRG is different from that of other local corre- 
lation methods (such as local CCSD) since these 
require the correction to the mean-field reference 
to be small and to possess finite correlation length. 

The DMRG ansatz possesses a further technical 
advantage. Since no localisable Hartree-Fock refer- 
ence is required, the favorable scaling of the DMRG 
is obtained in any local basis, and does not in 
particular need separate localisation in the occu- 
pied and virtual spaces. This is particularly ad- 
vantageous when modelling correlated states which 
possess a shorter quantum correlation length than 
their parent mean-field reference (e.g. systems with 
small Hartree-Fock bandgaps) , for which orbital lo- 
calisation may be more difficult. 

5. Long molecules: The Matrix Product State embod- 
ies a finite correlation length as measured along 
the DMRG orbital lattice, rather than as neces- 
sarily exists in the physical space. Consequently, 
we can only hold M constant as the system size 
increases - and obtain the same relative accuracy 
- if the locality in the physical system maps geo- 
metrically onto a one-dimensional lattice, i.e. the 
system is extended in only one of its three dimen- 
sions as in a long molecule is avoided. If this is not 
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the case, we will require M to scale exponentially 
in the width of the system to maintain accuracy, 
much like Full Configuration Interaction (FCI) or 
CASSCF theory. In practice, we have shown that 
with reasonable M we can still obtain highly accu- 
rate DMRG energies even in non-one-dimensional 
molecular systems with up to O(40) active orbitals 
- i.e. too large to treat using FCI theory - but to 
model much larger extended non-one-dimensional 
networks of strongly interacting electrons, further 
progress in the DMRG method will be required. 



III. A QUADRATIC SCALING PARALLELISED 
DMRG ALGORITHM 



The full computational scaling of a single conventional 
DMRG sweep is 0(M 2 k 4 ) + 0(M 3 k 3 ). Here, the 0(k 4 ) 
scaling arises in essence from the number of two electron 
integrals Vijki in the Hamiltonian H, written in second 
quantisation as 



H = ^2 '</",".. + ^^Vijuala^akdi. (HI.l) 

ij ijkl 



Recall that M can be kept fixed, independent of system 
size in a long molecule. Thus to implement a quadratic 
scaling DMRG algorithm we need only screen the con- 
tributions from the 2 electron integrals. This can be 
achieved by working in a localized basis. (Note we can 
use any localised orthonormal basis and we do not need 
to separately localise the occupied and virtual spaces as 
is commonly required in local correlation methods. For 
example, later in this work, we shall use the basis of over- 
lap symmetrically orthonormalised atomic orbitals). As 
is well understood, in a large system described in a lo- 
calised basis, the number of significant two-electron inte- 
grals below a given threshold scales only quadratically as 
non-classical Coulomb integrals i.e. integrals of the form 
Vijkl = §(i(l)Z(l)|j(2)fc(2)) where i(l),i(l) or j(2),k(2) 
functions are widely separated, vanish exponentially with 
the separation between i, I or j, k centres. 

In the DMRG, we work with a number of intermedi- 
ate combinations of operators on each of the blocks of 
orbitals which are subsequently combined to construct 
the Ml H IHliH. A DMRG sweep, consisting of Oik) 
sweep iterations (each comprising a different block con- 
figuration), requires 0(M 2 k 4 )+0(M 3 k 3 ) time, 0(M 2 k 2 ) 
memory, and 0(M 2 k 3 ) disk storage. These asymptotic 
costs originate from manipulating the two-index interme- 



FIG. 1: Standard block configuration in DMRG. From left to 
right, L, o L , o R , R. 
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diate operators Ay , Py , Py , Qy on the various blocks: 



.blk _ 


CLiCLj 


(III.2) 


nblk 

^•(Gblk) - 




(IH.3) 


pblk _ 

Mj(gblk) - 


y^ Vijuako-i 

fcZGblk 


(IH.4) 


^y(^blk) — 


^ Xi kia k a i 

feieblk 


(III.5) 


Xijkl — 


Vijkl — Vjiki — Vijik + Vjuk- 


(III.6) 



To begin, we employ screening to determine a set of 
significant two index operators that must be considered 
on each block, according to the following criterion: 

^j^gbik) discard if v^ki < threshi for all kl £ blk 

^■yfebik) discard if Xijki < thresh! for all kl £ blk 

PijOt blk) discard if Uyjy < threshi for all kl G blk 

Q^(rfi,ik) discard if x^ki < threshi for a ll kl 6 blk. 

(HL7) 

In a DMRG block configuration, there are four kinds of 
blocks: the left block L, an orbital to be blocked with the 
left block o £ , an orbital to be blocked with the right block 
0^, and the right block R (Fig. QJ. Without screening, 
the number of two-index operators that must be con- 
sidered on each block is 0(k 2 ), but in each case this is 
reduced to 0(k) after screening since eqns. ()IIL7() require 
centres i,j to be close in space. Since the number of op- 
erators is reduced, we also reduce the memory cost to 
0(M 2 k) per sweep iteration (block configuration). The 
disk usage is reduced to 0(M 2 k 2 ) per sweep. 

Next, we consider the computational costs of the dif- 
ferent manipulations involving the two index operators in 
each of the three stages of a sweep iteration: (1) blocking, 
(2), solving for the wavefunction, and (3) decimation: 

1. Blocking: Here we construct representations of the 
operators in the tensor product space of a large 
block and an additional orbital; for concreteness, 
we take the large block as the left block L, and 
the additional orbital as o^, and we consider the 
operator Py. First, we accumulate P^j,P° L in the 
new space {L} ® {°l}; for each such term, the ac- 
cumulation requires 0(M 2 ) time. Since there are 
0(k) screened Py operators on both blocks L and 
o L , in total this requires 0(M 2 k) time per blocking 
step and thus 0(M 2 k 2 ) time per sweep. Next, we 
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sum over the new terms appearing in eqns. (|III.2|I . 
(|111.3|l that arise from the combinations ft G L, I 6 
o L and ft G o L .l g L; each such term requires 
0(M 2 ) time per blocking step. Without screen- 
ing the number of new terms become O(ft), but 
with screening we discard any contributions where 
Vijki < thresh2 for all kl G blk and this decreases 
the number of new terms to 0(1) for each signifi- 
cant Pij operator per blocking step. Consequently, 
the time to accumulate the additional contributions 
is 0{M 2 ) x O(l) x no. significant P tj = 0(M 2 k) 
per blocking step, or 0(M 2 k 2 ) time per sweep. Re- 
peating this analysis for the j4y, fly, Qij operators, 
we observe that these also involve 0(M 2 k 2 ) time 
per sweep. 

2. Solving for the wavefunction: In an iterative David- 
son algorithm 0] , the contributions of , Qij to 
the Hamiltonian matrix multiply takes the form 

Each CsS requires 0(M 3 ) time, and thus the over- 
all cost is determined by the number of ij indices 
to sum over. From the screening criterion threshi, 
this is 0(k) for each block configuration, and thus 
the total time for a single Hamiltonian multiply 
takes 0(M 3 k), or 0(M 3 k 2 ) per sweep. 

3. Decimation: In the decimation for each two index 
operator, each transformation takes 0(M 3 ) time. 
After screening, only 0(k) ij indices need be con- 
sidered per block, and thus the time to transform 
all Aij, Bij, Pij, Qij operators is 0(M 3 k) per renor- 
malisation step, or 0(M 3 k 2 ) per sweep. 



FIG. 2: Geometries of chemical systems used in the LDMRG 
study in sec. II VI 
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performance of the quadratic-scaling algorithm and ro- 
bustness of the screening criteria, (iii) convergence of the 
LDMRG ansatz, and (iv) errors compared against stan- 
dard correlation methods. We have chosen two classes 
of systems as representative "long" molecules: planar 
all-trans-polyenes CkHk+2 ranging from ft = 4, 8, ... ,48 
(modelled in the ^-active space) and hydrogen molecule 
chains (Ha)fc/2 ranging from ft = 10, 20, . . . , 100. The ge- 
ometries of the polyenes (based on 47]) and hydrogen 
chains are given in fig. We note that although the 
bond-lengths in the hydrogen molecule chains are alter- 
nating, the molecules are still spaced sufficiently closely 
to be interacting. 



In summary, integral screening in the LDMRG reduces 
the total computation cost per sweep to 0(M 3 k 2 ) + 
0(M 2 k 2 ) time (i.e. quadratic scaling, since M in long 
molecules is independent of system size for a chosen ac- 
curacy and hence a constant), 0(M 2 k) memory, and 
0(M 2 k 2 ) disk. Table [I] summarises the key operations 
and costs of the screened algorithm. 

Finally, we note that the above screening procedure 
is easily combined with the parallelised algorithm em- 
ployed in our previous calculations [3 If . Once the list of 
screened ij indices is determined via eqns. IIII.7I the sig- 
nificant operators are distributed over the processors, and 
all manipulations involving these operators are then car- 
ried out in parallel. This screened parallelised algorithm 
has been employed to perform the calculations described 
in the current work. 



IV. NUMERICAL ANALYSIS OF THE LDMRG 
IN LONG SYSTEMS 

In the current section, we report our numerical in- 
vestigations of (i) the accuracy and extensivity of the 
LDMRG ansatz in long molecules, (ii) computational 



A. Computational details 

All electronic integrals were obtained using the Psi3.2 
package [isj. We used an STO-3G minimal basis for the 
polyene calculations and an STO-6G minimal basis for 
the hydrogen molecular chains [i^. l5fj| . Polyene calcu- 
lations were performed in the 7r z -active space spanned 
by one p z orbital on each carbon center, with each car- 
bon atom contributing one electron; thus in CkHk+2 
we used a (ft, ft) active space. The remaining electrons 
were placed in doubly occupied Restricted Hartree-Fock 
(RHF) orbitals generated by the PSI3.2 program. Cal- 
culations on the hydrogen chains correlated all electrons. 

We used a localized orthonormal basis as input to the 
LDMRG calculations. In both the polyenes and hydro- 
gen chains this was obtained by symmetrically orthonor- 
malising (S^ 1 / 2 ) the atomic orbital basis. The orthonor- 
malised orbitals were then ordered in their natural topo- 
logical order i.e. in the order of their originating atoms 
along the chain. Since each atom contributes only one 
basis function, this ordering is unique. 

The LDMRG calculations were performed with the 
parallel Block code [3l| with integral screening as de- 
scribed in sec. IIIII on 4-18 processors. Except where 
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TABLE I: Time, memory, and disk costs associated with the two-index operators in the original DMRG and screened LDMRG 
algorithms. The two-index operators determine the asymptotic computational costs of the algorithm. 



Operator 


Blocking 


Solving 


Decimation 


Memory 


Disk 




DMRG LDMRG 


DMRG LDMRG 


DMRG LDMRG 


DMRG LDMRG 


DMRG LDMRG 


A - R 


M 2 k 3 M 2 k 2 


M 3 k 3 M 2 k 2 


M 3 k 3 M 3 k 2 


M 2 k 2 M 2 k 


M 2 k 3 M 2 k 2 


Pij i Qij 


M 2 k 4 M 2 k 2 


M 3 k 3 M 2 k 2 


M 3 k 3 M 3 k 2 


M 2 k 2 M 2 k 


M 2 k 3 M 2 k 2 



stated otherwise (see sec. IIV C|) . we applied screening 
thresholds of threshi = l(T 7 E h and thresh 2 = KT 20 E h . 
No spatial symmetry was used. DMRG sweeps were per- 
formed with progressively increasing M values (a sweep 
schedule) and a small amount of random noise (between 
1CP 6 — 1CP 9 in the matrix norm) was added to the density 
matrix in the early sweeps (M ^ 100) to prevent loss of 
quantum numbers |29l l31|. A typical schedule to obtain 
AI = 50, 100, 250 DMRG energies is as follows: sweeps 
1-6: M = 50 (with noise), sweeps 7-12: M — 50, sweeps 
13-18: M = 100 (with noise), sweeps 19-24: M = 100, 
sweeps 25-30: M = 250. We have converged our LDMRG 
energies to 8 significant figures; unconverged digits are 
denoted in italics. Because of the complexity of the ab 
initio DMRG method and the non-linearity of the opti- 
misation, there is a small dependence of the DMRG en- 
ergies on the precise computational setup (e.g. the way 
in which M is increased in sweeps) which may lead to 
some variation in the last significant digit. 

B. Accuracy and Extensivity of the DMRG ansatz 



FIG. 3: All-trans-polyenes: Active space correlation energy 
from MP2, CCSD, CCSD(T) and LDMRG(250) as a func- 
tion of polyene chain length. On the scale of the graph, the 
LDMRG and CC results nearly overlap. 
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In Tables HT1 and ITTT1 we present the energies obtained 
with our quadratic scaling LDMRG algorithm for the 
all-trans-polyene series and hydrogen molecular chains. 
For comparison, we also present second-order M0ller- 
Plesset (MP2) and Coupled Cluster calculations (CCSD, 
CCSD(T)) obtained using the Psi3.2 (hydrogen chains) 
and Dalton 2.0 [5l| (active-space polyenes) packages. 

In the largest M LDMRG calculations shown, the cor- 
relation energies are exact correlation energies for the 
many-particle Schrodinger equation to the digits dis- 
played. For example, in the polyenes, calculations at the 
LDMRG(500) and LDMRG(IOOO) level did not change 
the energy in the /zEh-range. To confirm the exactness 
of our LDMRG calculations, we also performed explicit 
active space FCI calculations (using Molpro 2002.6 
[H2|) f° r C4H 6 and CeH 8 and obtained agreement to all 
displayed digits. Following the discussion in sec. 

MM 

the hydrogen molecular chain energies are presented to 
10/iEh precision, corresponding to eight significant fig- 
ures in the electronic energy of the longer chains. There 
are only improvements of the order of 1/xEh when going 
to LDMRG(IOO) and thus the LDMRG(50) correlation 
energies for the hydrogen molecular chains are exact to 
the digits displayed. 

The largest Hilbert space considered (for the (H 2 )5o 
system containing 100 electrons in 100 orbitals) has 



dim(7Y) = 10 58 . That we are able to obtain a numerically 
exact correlation energy with the LDMRG illustrates the 
compactness of the LDMRG description in systems that 
are still interacting but have finite correlation lengths, 
which allows us to keep M fixed as the system size grows 
(see sec. Ill A|) . A related feature of the LDMRG ansatz is 
that of size-consistency/extensivity of the energy, which 
we now discuss. 

In fig. |21 we plot the active space correlation energy 
Ecorr.act as a function of polyene chain length. A clear 
linear relationship between chain length and correlation 
energy is observed. Fig. 0] shows in detail how the ac- 
tive space correlation energy as well as the total energy 
Etot per additionally introduced C4H4-unit converges to 
a constant in the limit of long polyenes. 

We also performed a series of lower accuracy LDMRG 
calculations for the polyenes, with M = 5 — 40 states. 
Due to the variational nature of the DMRG, these ap- 
proach the exact energy from above. Fig. presents the 
logarithm of the percentage error in the correlation en- 
ergy relative to the "exact" LDMRG(250) energies, as a 
function of chain length. In small systems the LDMRG 
calculations are exact, since the DMRG states span the 
whole A^-particle space. In the longer polyenes, the per- 
centage errors increase to a saturating value, demonstrat- 
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FIG. 4: All-trans-polyenes: Exact (LDMRG(250)) total en- 
ergy per additionally introduced CiHVunit. The inlay shows 
the active space correlation energy at CCSD, CCSD(T) and 
LDMRG(250) level of theory per additionally introduced 
C4H4-unit. 
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FIG. 5: All-trans-polyenes: Relative errors in the active space 
correlation energies for LDMRG with various M (compared 
to the exact LDMRG(250) results). The black marked curves 
are the errors of MP2, CCSD, and CCSD(T) as reference. The 
plot shows absolute magnitudes in logarithmic scale. 
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ing the size-extensivity of the approximate LDMRG cal- 
culations. Similar observations can be made for the hy- 
drogen molecular chains. 



C. Computational Scaling and Screening 
Robustness 

In Figure El we present the asymptotic computational 
scaling of the sweep time for the LDMRG calculations 
as a function of the number of active orbitals for the 



polyenes and hydrogen molecular chains. Here sweep 
times were measured after several sweeps at a given M 
level had been performed, to remove the bias that occurs 
immediately after a transition from a lower M calculation 
in the sweep schedule. 

We fitted the timing data to obtain the computational 
scaling of LDMRG as a function of the number of active 
orbitals. The scaling exponents for the polyenes and the 
(H2)fe/2-chains with different M and screening thresholds 
are given in tab. [3 

In the polyenes we find a reduced scaling of near- 
quadratic order, with an exponent between 2.1 — 
2.2. For reasonable screening thresholds (i.e. threshi 
10 _6 Eh— 10 _8 Eh) no significant differences in the scaling 
is observed. We also do not see a significant scaling de- 
pendence on M. In the hydrogen chains a similar reduced 
scaling was found, in this case with exponents ranging 
from 2.2 — 2.4. In both cases, it is clear that the screened 
LDMRG algorithm has reduced the computational scal- 
ing of the DMRG to quadratic order. As an example of 
absolute times per sweep, for the largest system (H 2 )5o 
using 18 2.0 GHz Opteron processors, we required 27min 
for M = 50, 37min for M = 100, and 73min for M = 250. 

Since the LDMRG employs screening, we should assess 
the robustness of the criterion that is used. To this end, 
we studied the polyene correlation energies computed 
with screening thresholds (threshi) of 10~ 6 E h , 10- 7 E h , 
10 _8 Eh and 10 _20 Eh (the energy of the latter can be con- 
sidered unscreened). A selection of results is presented 
in tab. IIVI We observe the correlation energy to be 
converged at the /zEh level with respect to the screening 
threshold when threshi = 10~ 7 Eh, which is the reason 
for using this setting during this study. In practice, this 
threshold could be relaxed for lower accuracy calcula- 
tions. 



D. Sweep and Error Convergence 

There are two types of convergence in DMRG calcu- 
lations. The first is the convergence of the energy as a 
function of the number of sweeps, holding the number 
of DMRG states M fixed. We observed that on average, 
convergence was achieved in only 4—6 sweeps for small M 
and 2 — 4 sweeps for large M values (not inclusive of the 
noise sweeps and the preceding sweeps in the schedule). 

The second type of convergence relates to the approach 
of the DMRG energy to the exact energy as the number of 
retained states M is increased. Here, we analyse our data 
for the DMRG calculations on polyenes with different M 
values. The precise analytic form of the DMRG energy 
convergence as a function of M has been a matter of 
debate in the literature |H |H |M El HI- We have 
previously found good agreement with the proposed form 
of Okunishi et al. |53J, |5J] , 



AE(M) - exp [-k (logM) c 



a = 2 



(iv.i) 



which is slower than exponential but still faster than al- 
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TABLE II: All-trans-polyenes: Dimension of the FCI determinant space, total RHF energy, RHF active space electronic energy; 
active space correlation energies at MP2, CCSD, CCSD(T) and different LDMRG levels of theory. All energies are given in 
hartrees. 



Molecule dim(H) a 


^RHF 




ERHF,el,act 


Ecorr, 


act 


























MP2 


CCSD 


CCSD(T) LDMRG(50) LDMRG(IOO) LDMRG(250) C 


C4H5 


3 gxlO 1 


-153.006 


364 


_3 


.169 


490 


-0. 


.046 


529 


-0 


.091 


435 


-0.091 


668 


-0.091 


502 


COTIV ^ 




COTIV 




Cs H10 


4.9 xlO 3 


oc\ a on 

-304.889 


389 


-8 


.426 


391 


-0. 


.090 


346 


-0 


.176 


445 


-0.177 


797 


-0.177 


127 


COTIV 




COTIV. 




Ci -Hi „ 
W2-tll4 


s ^ v 1 n 5 
0.0 X 1U 


-456.773 


412 


1 1 


.ooy 


000 


-0. 






-u 


9 fin 


77Q 


n 9fi3 


1 


-U.ZDZ 




f) 9fi9 
-U.ZDZ 


9Q7 


COTIV. 




CieHis 


1.7x10 s 


-608.657 


556 


-21 


.345 


452 


n 
-u 


.178 


366 


-0 


.345 


003 


-0.349 


327 


-0.347 


399 


-0.347 


403 


COTIV. 




C20H22 


3.4xl0 10 


-760.541 


718 


-28 


.542 


181 


-0. 


.222 


434 


-0 


.429 


210 


-0.435 


082 


-0.432 


490 


-0.432 


498 


COTIV. 




C24H26 


7.3xl0 12 


-912.425 


883 


-36 


.090 


721 


-0 


.266 


507 


-0 


.513 


414 


-0.520 


840 


-0.517 


579 


-0.517 


591 


COTIV. 




C28H30 


1.6xl0 15 


-1064.310 


048 


-43 


.931 


953 


-0 


.310 


582 


-0 


.597 


618 


-0.606 


599 


-0.602 


668 


-0.602 


684 


conv. 




C32H34 


3.6xl0 17 


-1216.194 


214 


-52 


.023 


816 


-0. 


.354 


658 


-0 


.681 


822 


-0.692 


358 


-0.687 


757 


-0.687 


777 


COTIV. 




C36H38 


8.2xl0 19 


-1368.078 


379 


-60 


.334 


842 


-0. 


.398 


734 


-0 


.766 


027 


-0.778 


118 


-0.772 


846 


-0.772 


870 


COTIV. 




C40H42 


1.9xl0 22 


-1519.962 


544 


-68 


.840 


593 


-0 


.442 


810 


-0 


.850 


231 


-0.863 


879 


-0.857 


935 


-0.857 


962 


-0.857 


963 


C44H46 


4.4xl0 24 


-1671.846 


710 


-77 


.521 


543 




e 






□ 








-0.943 


024 


-0.943 


055 


-0.943 


056 


C48H50 


l.OxlO 27 


-1823.730 


875 


-86 


.361 


727 




□ 










□ 




-1.028 


113 


-1.028 


147 


-1.028 


149 



\fkf3 
.1 \N 



Nft' 



here dim(W)= 



( k/2y 

\k/2) 



No point group 



a dim(7i)= 
symmetry used. 

'The active space electronic energy contains the core-active 
Coulomb and exchange contributions, but no nuclear repulsion. 
C AU calculations with M 250 converged. 

d " conv." denotes converged results, where increased M did not 
change the significant figures in the energy. 

e C44H4g and C48H50 could not be computed by Dalton due to 
an address limitation. 



TABLE III: (H2)fc/2-chains: Dimension of the FCI determinant space, total RHF energy, RHF electronic energy; correlation 
energies at MP2, CCSD, CCSD(T) and LDMRG(50) levels of theory. All energies are given in hartrees. 



E CO rr 



Molecule dim(H) a 


Erhp 


ErhF.cI 




MP2 


CCSD 


CCSD(T) LDMRG(50)'' 


(Ha) 8 


6.4xl0 4 


-5.553 


26 


-16.036 48 


-0. 


.068 


34 


-0.101 


93 


-0.102 


04 


-0.102 09 


(h 2 ; 


)io 


3.4xl0 10 


-11.088 


22 


-38.784 11 


-0. 


.137 


53 


-0.203 


77 


-0.204 


04 


-0.204 15 


(h 2 ; 


)15 


2.4xl0 16 


-16.623 


18 


-64.210 83 


-0. 


.206 


72 


-0.305 


61 


-0.306 


03 


-0.306 21 


(h 2 ; 


>20 


1.9xl0 22 


-22.158 


14 


-91.378 72 


-0. 


.275 


91 


-0.407 


44 


-0.408 


02 


-0.408 26 


(h 2 ; 


>25 


1.6xl0 28 


-27.693 


11 


-119.841 65 


-0. 


.345 


10 


-0.509 


28 


-0.510 


01 


-0.510 32 


(h 2 ; 


)30 


1.4xl0 34 


-33.228 


07 


-149.336 70 


-0. 


.414 


29 


-0.611 


12 


-0.612 


01 


-0.612 38 


(h 2 ; 


)35 


1.3xl0 40 


-38.763 


03 


-179.690 11 


-0 


.483 


48 


-0.712 


95 


-0.714 


00 


-0.714 44 c 


(h 2 ; 


)40 


1.2xl0 46 


-44.297 


99 


-210.778 35 


-0. 


.552 


67 


-0.814 


79 


-0.815 


99 


-0.816 49 


(h 2 ; 


)45 


l.lxlO 52 


-49.832 


95 


-242.509 08 


-0. 


.621 


87 


-0.916 


63 


-0.917 


98 


-0.918 55 


(h 2 ; 


)50 


l.OxlO 58 


-55.367 


92 


-274.810 58 


-0. 


.691 


06 


-1.018 


47 


-1.019 


98 


-1.020 61 



a No point group symmetry used. 

'All calculations with M ^ 50 converged. 

c This result lies very close on the rounding border, LDMRG(50) 
rounds the last digit to 4, LDMRG(IOO) to 5. 



gebraic. In fig. [7|we plot the logarithm of the percentage 
error in the correlation energy against (logM) 2 , which 
shows a clear linear fit. By contrast, the inlay plots 
the logarithm of the percentage error against M, which 
demonstrates that the error indeed does not decay ex- 
ponentially. Fitting our data (omitting M = 5) to the 
general form of eqn. IIV.1I we obtained an exponent of 
a ~ 1.6 — 1.8. Fixing a = 2, we obtain values between 



k = 1.80 ± 0.03 (for C12H14) and k = 1.45 ± 0.03 (for 
C48H50). 

Corresponding to the rapid energy convergence we also 
observed a rapid decrease of the truncated weight of the 
density matrix as M is increased. This shows that the 
local representation is well suited to the chemical system 
and physical problem at hand [35L I45I l55| . 
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TABLE IV: All-trans-polyenes: Active space correlation en- 
ergies from LDMRG(250) with screening thresholds 10" 6 E h 
and 10 _7 Eh; absolute and relative errors of 10" 6 E h-screemng 
(compared to the exact results from 10 _7 Eh-screening). All 
energies are given in hartrees. 



E CO rr,act 



Molecule 


(10~ 6 E h ) 


(10- 7 E h ) 


A abs 


A rcl [%] 




CgHio 


-0.177 127 


conv. 










CieHig 


-0.347 405 


-0.347 404 


-0.000 001 


-0.46x10 


-5 


C32H34 


-0.687 765 


-0.687 777 


0.000 008 


1.52x10" 


-5 


C40H42 


-0.857 942 


-0.857 963 


0.000 021 


3.0J xl0~ 


-5 


C48H50 


-1.028 093 


-1.028 149 


0.000 056 


6.4J xl0~ 


-5 



FIG. 6: (H2)fc/2-chains (circles) and all-trans-polyenes 
(squares): Asymptotic timing data (i.e. total time per sweep) 
of LDMRG with M = 50 (filled marks) and M = 250 (un- 
filled marks) for 10 _7 Eh-screening in log-log-representation 
with linear fit. 
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E. Comparison with Perturbation and Coupled 
Cluster Theories 

Since with our larger M LDMRG calculations we have 
exact energies at our disposal we can analyze the errors 
at the various levels of theory in more detail. 

In the polyene calculations (tab. [HJ the largest DMRG 
absolute error is 35/iEh for the C48H50 molecule at the 
AI = 50 level. This corresponds to ~ 10~ 3 % of the ex- 
act active space correlation energy, and ~ 1CP 5 % of the 
exact total active space electronic energy. Compared to 
the Coupled Cluster errors, LDMRG(50) is already bet- 
ter by 2-3 orders of magnitude. The LDMRG(IOO) gives 
a further order of magnitude improvement, and is essen- 
tially exact. In our more approximate calculations (fig. 
13 we find that LDMRG with M = 10 performs better 
than MP2, and with M = 15 better than CCSD and 
CCSD(T). The results for M = 5 are not reliable due to 
loss of important quantum numbers. 



FIG. 7: All-trans-polyenes: Convergence of the relative er- 
rors in the active space correlation energies for LDMRG as a 
function of M (compared to the exact LDMRG(250) results). 
The main plot shows magnitudes in logarithmic scale over 
log (M) 2 (with linear fit), the inlay shows them over log(M). 




Ho ' ls ' 2jo ' i!i ' 3X1 ' i!i ' Zo 

log(iW) 2 



Surprisingly we observe that the CCSD(T) results lie 
below the exact energies computed with LDMRG. Al- 
though CCSD(T) is not variational in general, it is still 
uncommon to obtain an energy below the exact energy at 
an equilibrium geometry. In general the triples correction 
performed relatively badly for the polyenes. 

In case of the hydrogen chains, the convergence of the 
LDMRG with M was more rapid and results were already 
exact with M = 50. CCSD(T) also performed better in 
this system, the triples correction improved on CCSD by 
1/2 order of magnitude, and the resulting energies were 
consistently above the exact energies. 

V. THE METAL-INSULATOR TRANSITION IN 
LINEAR HYDROGEN 

As an example of a challenging electronic problem, we 
studied the symmetric and asymmetric bond stretching 
in a linear Hso-chain. In both these cases, the system 
transitions from a state with metallic correlations at com- 
pressed geometries to an insulating state with strong mul- 
tireference correlation in the dissociation region. This 
bond breaking process hence exhibits a varying nature of 
chemical bonding and electron correlation. 

In case of the symmetric dissociation we begin with a 
uniform bond-distance between all H-atoms of R=1.0ao, 
and stretch all 49 bonds symmetrically and simultane- 
ously to R=1.2, 1.4, 4.2ao- The final structure con- 
sists of 50 equidistant, nearly-independent H-atoms on a 
line. 

In case of the asymmetric dissociation we distinguish 
alternating bonds as intermolecular and intramolecular 
with Ri n ter and Rmtra- The first geometry is Ri n tra=Rintor 
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TABLE V: Asymptotic scaling exponents (with standard error) of LDMRG depending on M and the screening threshold. 



Scaling exponent 



System (threshi [E h ]) LDMRG(50) LDMRG(IOO) LDMRG(250) 
C fc Hfc + 2 (1(T 6 ) 2.12±0.02 2.11±0.01 2.07±0.09 

C fc H fc+ 2 (10 -7 ) 2.11±0.02 2.11±0.01 2.10±0.03 

CfcH fe+ 2 (1(T 8 ) 2.12±0.01 2.07±0.02 2.09±0.03 

C fe H fc+ 2 (KT 20 ) 3.27±0.08 3.33±0.10 3.53±0.06 

(H 2 ) fc/ 2 (1CT 7 ) 2.36+0.06 2.18±0.05 2.16±0.04 



—IAslq. In the following geometries Rintra is kept fixed 
at IA&q while Rinter grows to Ri ntor =1.6, 1.8, 4.2ao- 
The final structure consists of 25 equidistant, nearly- 
independent H2-molecules at equilibrium bond distance 
on a line. 

We computed the electronic energy using the LDMRG 
(with up to 1000 states) in the minimal STO-6G basis, 
where we correlated all 50 electrons in 50 orbitals. 

All calculations were carried out in the STO-6G basis 
correlating all electrons (50 electrons in 50 orbitals). The 
LDMRG calculations again used the S^ 1 ^ 2 basis. 



A. Symmetric Dissociation 

The calculated energies for the symmetric dissociation 
are summarized in tab. IVII The potential energy curves 
at RHF, MP2, and exact LDMRG level of theory are 
plotted in fig. [S] It can immediately be seen how the 
contribution of correlation increases along the dissocia- 
tion coordinate: In the dissociation limit the share of the 
correlation energy in the total energy grows to ~20% and 
in the electronic energy to ~7%, which emphasizes the 
importance of nondynamic correlation in this problem. 

As is expected, RHF and MP2 behave poorly as the 
chain dissociates. The Coupled Cluster energies cannot 
even be converged for bond lengths R>2.0ao. This is 
a fundamental problem in CC theory that is well docu- 
mented e.g . in the work of Takahashi and Paldus and oth- 
ers |56l l57l IHs[ where in one-dimensional systems, even 
for physically relevant coupling parameters, the Coupled 
Cluster doubles equations may have no real solutions. 
The correlation energy errors for different methods rela- 
tive to the exact LDMRG results are shown in Fig. |5| 

It is understood that we need to retain more states in 
the LDMRG in the metallic regime if we start from a local 
atomic orbital basis, since we need to capture the delocal- 
isation and long-range off-diagonal correlations ^5|. We 
find that both the convergence with the number of sweeps 
as well as with M is slower as compared to calculations 
in the nonmetallic regime. At R=1.0ao LDMRG(50) 
is worse than CCSD, LDMRG(IOO) slightly worse than 
CCSD(T), and for R<1.6a LDMRG(50) is still worse 
than CCSD(T). In the metallic region LDMRG required 
M = 500 to converge to the numerically exact result. 
In essence, by using orthonormalised atomic orbitals, we 



FIG. 8: Symmetric dissociation of H50: Potential energy 
curves from RHF, MP2, and exact LDMRG. On the scale of 
the graph, the few available CCSD and CCSD(T) datapoints 
were indistinguishable from the LDMRG data. 



♦ — RHF 
■ — MP2 

• — LDMRG(500) 




R[aJ 



FIG. 9: Symmetric dissociation of H50: Relative errors in the 
correlation energies at MP2, CCSD, CCSD(T), and different 
LDMRG levels of theory (compared to the exact LDMRG 
results) in logarithmic scale. 
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TABLE VI: Symmetric dissociation of H50: Total RHF energy, RHF electronic energy; correlation energies at MP2, CCSD, 
CCSD(T) and various LDMRG levels of theory. All energies are given in hartrees. 







Ecorr 


R [ao] Erhf 


ErhP,c1 


MP2 CCSD CCSD(T) LDMRG(50) LDMRG(IOO) LDMRG(250) LDMRG(500) a 



1 





-16 


864 


88 


-191 


825 


14 


-0 


361 


45 


-0.407 


29 


-0 


417 


39 


-0.402 72 


-0 


417 


27 


-0.419 14 


-0.419 19 


1 


2 


-22 


461 


27 


-168 


261 


49 


-0 


401 


83 


-0.470 


11 


-0 


483 


30 


-0.475 90 


-0 


485 


21 


-0.486 35 


-0.486 38 


1 


4 


-25 


029 


76 


-150 


001 


38 


-0 


444 


73 


-0.543 


03 


-0 


559 


36 


-0.557 16 


-0 


563 


30 


-0.564 00 


-0.564 02 


1 


6 


-26 


062 


25 


-135 


412 


42 


-0 


491 


88 


-0.631 


18 


-0 


650 


89 


-0.652 72 


-0 


656 


74 


-0.657 18 


-0.657 19 


1 


8 


-26 


265 


98 


-123 


466 


13 


-0 


545 


50 


-0.741 


67 


-0 


765 


47 


-0.769 82 


-0 


772 


42 


-0.772 66 


-0.772 67 


2 





-26 


008 


20 


-113 


488 


34 


-0 


607 


89 


-0.883 


29 


-0 


912 


70 


-0.916 11 


-0 


917 


76 


-0.917 89 


conv. 


2 


4 


-24 


835 


76 


-97 


735 


87 


-0 


768 


83 




C 






□ 


-1.324 16 


-1 


324 


77 


-1.324 81 


conv. 


2 


8 


-23 


360 


81 


-85 


846 


62 


-0 


995 


30 




□ 






□ 


-1.913 81 


-1 


913 


98 


-1.913 99 


conv. 


3 


2 


-21 


896 


33 


-76 


571 


41 


-1 


307 


78 




□ 






□ 


-2.671 90 


-2 


671 


95 


conv. 


conv. 


3 


6 


-20 


574 


29 


-69 


174 


36 


-1 


723 


32 




□ 






EH 


-3.528 46 


-3 


528 


48 


conv. 


conv. 


4 


2 


-18 


955 


95 


-60 


613 


15 


-2 


558 


99 










□ 


-4.793 76 




conv. 


conv. 


conv. 



a All calculations with M ^ 500 converged. 

conv." denotes converged results, where increased M did not 
change the significant figures in the energy. 

c The Coupled Cluster calculations could not be converged. See 
text. 



are starting from a particularly unfavourable one-particle 
basis to describe metallic behaviour. By performing the 
DMRG in a set of separately localised occupied and vir- 
tual orbitals such as Boys orbitals , we expect that the 
degradation in efficiency of the DMRG would be avoided. 



B. Asymmetric Dissociation 

The calculated energies for the asymmetric dissocia- 
tion are summarized in table IvTTl In this system, the re- 
stricted Hartree-Fock reference dissociates correctly (to a 
set of non- interacting hydrogen molecules) , which can be 
understood by changing to a localised basis in the space 
of restricted occupied orbitals. For this reason, the re- 
stricted MP2 and CC theories are also qualitatively cor- 
rect and we see that their energies (fig. II 01) lie parallel 
to the exact LDMRG values along the dissociation curve. 
Unlike in the symmetric dissociation, the correlation en- 
ergy saturates rapidly to ~ 1.8% of the total energy as the 
bonds are stretched. Fig. 1111 shows how the percentage 
errors in the correlation energy for the different methods 
decrease along the dissociation coordinate. Again, we 
see reduced performance of the DMRG in the metallic 
regime due to the unsuitability of the underlying orbital 
basis, but still a systematic convergence with M. For 
large Rintor we observed very rapid convergence with M 
and number of sweeps, and in fact for Ri nte r=4.2ao the 
LDMRG energy was already exact after 4 noise sweeps 
with M = 50. In the limit of complete dissociation, the 
CCSD theory becomes exact for this system and this is 
confirmed by convergence to the LDMRG results. 

In order to demonstrate the metal-insulator transition 
more explicitly we computed the one-particle reduced 



FIG. 10: Asymmetric dissociation of H50: Potential energy 
curves from RHF, MP2, and exact LDMRG. On the scale 
of the graph, the LDMRG, CCSD, and CCSD(T) curves are 
indistinguishable . 



1 — 1 — 1 — 1 — 1 — r 




l!s' ' ' 2.0' ' ' '2.5' ' ' '3.0' ' ' '3.5' ' ' 40 
R inter KJ 



density matrix 7 during our LDMRG calculations. In fig. 
IT2*1 we have plotted the off-diagonal decay of the a one- 
particle density matrix from element 7 25 25 — > 725 50- I n 
the metallic regime (short Rintor) we see the long-ranged 
oscillations in the off-diagonal elements, while in the insu- 
lating regime (long R^tcr) the off-diagonal elements de- 
cay much more rapidly. A similar picture is obtained 
from the density matrix during symmetric dissociation. 
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TABLE VII: Asymmetric dissociation of H50: Total RHF energy, RHF electronic energy; correlation energies at MP2, CCSD, 
CCSD(T) and various LDMRG levels of theory. All energies are given in hartrees. 

E 

-■-^corr 

Rintcr [a ] Erhf E RH F,oi MP 2 CCSD CCSD(T) LDMRG(50) LDMRG(IOO) LDMRG(250) LDMRG(500) a 



1.4 


-25.029 


76 


-150.001 


38 


-0.444 


73 


-0.543 


03 


-0.559 


36 


-0.557 16 


-0.563 30 


-0.564 00 


-0.564 02 


1.6 


-25.963 


71 


-142.811 


31 


-0.392 


61 


-0.516 


01 


-0.522 


30 


-0.523 02 


-0.523 64 


-0.523 67 


conv. b 


1.8 


-26.617 


68 


-136.573 


91 


-0.369 


20 


-0.505 


47 


-0.508 


73 


-0.509 38 


-0.509 48 


conv. 


conv. 


2.0 


-27.071 


82 


-131.089 


88 


-0.357 


01 


-0.503 


04 


-0.504 


97 


-0.505 48 


-0.505 50 


conv. 


conv. 


2.4 


-27.609 


24 


-121.878 


99 


-0.346 


20 


-0.507 


17 


-0.508 


03 


-0.508 37 


conv. 


conv. 


conv. 


2.8 


-27.873 


62 


-114.445 


23 


-0.341 


62 


-0.512 


77 


-0.513 


22 


-0.513 45 


conv. 


conv. 


conv. 


3.2 


-28.004 


68 


-108.324 


99 


-0.338 


67 


-0.516 


18 


-0.516 


42 


-0.516 56 


conv. 


conv. 


conv. 


3.6 


-28.069 


65 


-103.203 


54 


-0.336 


34 


-0.517 


50 


-0.517 


63 


-0.517 71 


conv. 


conv. 


conv. 


4.2 


-28.111 


00 


-96.924 


54 


-0.333 


73 


-0.517 


49 


-0.517 


54 


-0.517 58 


conv. 


conv. 


conv. 



"All calculations with M ^ 500 converged. 

b " conv." denotes converged results, where increased M did not 
change the significant figures in the energy. 



FIG. 11: Asymmetric dissociation of H50: Relative errors in 
the correlation energies at MP2, CCSD, CCSD(T), and differ- 
ent LDMRG levels of theory (compared to the exact LDMRG 
results) in logarithmic scale. 
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FIG. 12: Asymmetric dissociation of H50: Half-cross-section 
of the LDMRG a-one-particle density matrix at a-orbital no. 
25. 
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VI. CONCLUSIONS 

We began this work with the question of how to de- 
scribe nondynamic correlation in large systems with the 
restriction that such systems are large in only one di- 
mension. In our investigations, we have shown how the 
Density Matrix Renormalization Group (DMRG) pro- 
vides a natural answer to this problem. The Matrix 
Product State that underlies the DMRG is a local, varia- 
tional, size-consistent/size-extensive, and inherently mul- 
tireference ansatz that can efficiently exploit the spe- 
cial structure of quasi-one-dimensional correlation. Us- 
ing the intrinsic locality of the ansatz, we have formu- 
lated a quadratic scaling DMRG algorithm, using only 



a straightforward screening criterion without the impo- 
sition of correlation domains. With this active space 
method, we could then obtain numerically exact solutions 
of the many-particle Schrodingcr equation for all-trans- 
polyenes up to C48H50 (correlating the ^-electrons) and 
hydrogen molecular chains up to (112)50 (correlating 100 
electrons in 100 orbitals). 

By construction, a unique advantage of the LDMRG as 
compared to other local correlation methods is its ability 
to capture nondynamic correlation. We can take advan- 
tage of locality in multireference problems so long as the 
correlation length is finite. We have demonstrated the 
capability and efficiency of the LDMRG in these situa- 
tions by obtaining numerically exact correlation energies 
in the metal-to-insulator transition of linear Hso-chains, 
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where we correlate 50 electrons in 50 orbitals. 

With the possibility of accurately capturing nondy- 
namic correlation in long molecules, we can now begin to 
address the quantitative description of strongly interact- 
ing states as found in the spectrum of materials such as 
the conjugated organic polymers. Here, the natural next 
step would be to combine an LDMRG description of the 
nondynamic correlation in the active 7r-space with our re- 
cent developments in Canonical Transformation Theory 
|60|. to incorporate the dynamic correlation that arises 
in larger basis sets. 
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